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Abstract 


\ 

The  purpose  of  this  investigation  is  to  demonstrate  the  use  of 
boundary  element  techniques  for  the  dynamic  analyses  of  geometrical  ly 
repetitive  structures  using  the  traveling  wave  approach.  A  formulation 
of  the  boundary  element  method  (HEM)  for  2-D  isotropic  materials  is 
developed.  The  HEM  formulation  is  then  used  to  calculate  the  mass  and 
stiffness  matrices  of  one  bay  of  a  baseline  structure.  From  the  mass 
and  stiffness  matrices  a  transfer  matrix  is  developed  for  the  bay. 

Using  traveling  wave  theory,  the  transfer  matrix  is  then  used  to 
identify  the  dynamical  characteristics  of  a  multiple  bay  structure. 

Results  are  compared  to  continuum  theory,  f 
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DYNAMIC  STRUCTURAL  ANALYSES  USING  BOUNDARY  ELEMENT  METHODS 


I.  Introduction 


1.1  Background 

The  advent  of  high  speed  computers  has  led  to  a  revolution  in 
structural  analysis  methods.  Complex  numerical  algorithms  as  found  in 
finite  element  codes  are  new  routinely  used  in  the  analysis  of 
everything  from  coat  hangers  to  spacecraft.  As  structures  of  interest 
become  larger  and  more  complex,  more  accurate  and  efficient  algorithms 
are  developed.  Recently,  there  has  been  a  great  deal  of  interest  in  the 
use  of  boundary  element  methods  (BEM)  for  structural  analysis.  BEM 
offers  the  advantage  of  reducing  the  computational  size  of  the  problem 
compared  with  traditional  finite  element  methods  (FEM) .  The  application 
of  boundary  element  techniques  is  particularly  attractive  in  analyzing 
large  repetitive  truss-like  structures,  such  as  those  being  proposed  for 
orbiting  space  platforms.  A  typical  large  space  structure  will  be 
sensitive  to  wave  propagation  from  on-board  disturbances  such  as  gyros, 
actuators,  docking  procedures,  etc.  This  thesis  introduces  the  use  of 
boundary  element  theory  in  developing  the  wave  propagation  transfer 
matrix  for  two-dimensional  periodic  structures. 

1.2  Boundary  Element  Methods 

Boundary  element  theory  is  used  in  many  disciplines  and  has  proved  to 
be  an  efficient  and  elegant  method  for  solving  many  numerically 

I 
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intensive  problems.  The  boundary  element  method,  just  like  the  finite 
element  method,  is  based  on  the  approximate  solution  of  an  equation  or 
set  of  equations  describing  a  physical  problem.  Unlike  FEM  however,  BEM 
utilize  functions  that  identically  satisfy  the  governing  equations  and 
only  approximately  satisfy  the  boundary  conditions.  In  addition,  only 
the  boundary  of  the  given  problem  needs  to  be  discretized  when  using  the 
BEM.  This  greatly  reduces  the  modeling  effort  and  results  in  smaller 
matrices,  although  the  matrices  are  often  fully  populated.  A  further 
reduction  in  problem  size  is  accomplished  by  combining  a  traveling  wave 
approach  for  periodic  structures  with  BEM. 

1.3  Wave  Propagation 

The  study  of  wave  propagation  has  been  pursued  throughout  a  wide 
range  of  disciplines  including  solid  state  physics,  fluids,  power 
transmission,  etc.  Wave  propagation  theory  has  also  been  shown  to  be  a 
useful  tool  in  the  area  of  dynamic  structural  analysis.  Cremer  and 
Lielich  studied  flexural  motion  in  periodic  beams  (1) .  In  1964  Heckl 
defined  the  notion  of  propagation  coefficients  in  periodically 
supported,  undamped  grillages  (2:1335-1343).  More  recently,  von  Flotow 
employed  the  use  of  wave  propagation  theory  in  developing  a  transfer 
matrix  method  for  analyzing  periodic  structures  (3:509-519).  The 
transfer  matrix  method  requires  that  only  a  single  cell  of  the  truss 
structure  be  analyzed.  In  a  related  effort,  Signerolli  combined  the  use 
of  transfer  matrices  with  FEM  in  analyzing  a  two-dimensional  periodic 
truss  (4) . 
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1.4  Overview 


This  thesis  demonstrates  the  use  of  EEM  to  develop  the  transfer 
matrix  for  a  two-diirensional  periodic  truss.  Chapter  2  begins  with  the 
development  of  the  boundary  element  equation  for  a  two-dimensional 
isotropic  material.  The  boundary  element  equation  is  then  used  to 
develop  the  transfer  matrix  for  a  simple  periodic  structure  in  chapter 
3.  In  Chapter  4,  the  transfer  matrix  is  used  in  determining  the  wave 
propagation  behavior  of  a  baseline  beam  structure.  The  results  are  then 
compared  to  equivalent  continuum  models.  Chapter  5  discusses  the 
conclusions  and  recommendations. 
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II. 


2.1  The  Elastic  Dynamic  Equation 

The  governing  equation  for  the  behavior  of  a  deformable  and 
continuous  body  can  be  written  in  index  notation  as 


+  F,  -  nu,  (1) 

where  3  is  the  partial  derivative  operator,  a(j.  is  the  stress  tensor,  F,- 
is  the  body  force  vector,  m  is  the  mass  density,  and  is  acceleration. 
Another  useful  relationship  is  the  stress-strain  equation. 


a-.  =  2Ge--  +  \5--e  (2) 

ij  i  j  ljm  '  ' 

Where  G  is  the  shear  modulus,  eVj  is  the  strain  tensor,  <Sfj-  is  the 
Kronecker  delta,  and  X  is  Lame's  constant.  For  linear  isotropic  material 
properties,  strain  can  be  expressed  as 


=  h(  u,. 


i'i 


j  ’  > 


(3) 


where  Uj  is  the  displacement  vector.  Eqs.  (2)  and  (3)  can  be  combined 
to  give  stress  in  terms  of  displacement. 


!  aij  =  G(  *V j  +  uj'i  )  +  ^u'Vk 

i 

[ 


(4) 
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Eqs.  (1) ,  (2)  and  (3)  represent  15  different  equations  in  15 
unknowns;  6  a's,  6  e's  and  3  u's.  When  combined,  these  three  equations 
fully  describe  the  behavior  of  an  elastic  domain  and  will  be  used  to 
develop  the  boundary  element  equation  (5,210-211). 

To  formulate  the  solution  of  the  elasticity  equations,  Eq.  (1)  will 
be  multiplied  by  an  arbitrary  function,  u^,  and  integrated  throughout 
the  domain.  Body  forces  can  be  neglected. 

|  (  imii)dn  =  0  (5) 

This  equation  is  integrated  by  parts  until  all  of  the  differentials  are 
on  the  arbitrary  function,  u,.,  rather  than  the  unknown  displacements. 

Use  the  identity 


3xi 


Vij 


UjC* 


u  ,i 


(6) 


to  integrate  the  first  term  by  parts 


-  u .'fj..  -  Ujinu,.}  dfl  =0 


(7) 


Using  the  Guass-Divergence  theorem,  the  volume  integral  of  the  first 
term  may  be  written  as  a  surface  integral. 
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Jr  Wij)  ^ 


(8 


were  n-  is  the  direction  cosine  of  the  outward  unit  normal  vector.  Now 
substitute  Eq.  4  under  the  second  integral. 


Jrnj(Vij)  dr  -  £  +  uj-m)  +  +  u^}  dn  =  o  (9 


or 


..)  dr  - 


£  *  -  |n  °i.igui.i  * 

-  I  fli.M-k  an  -  f  an  =  o 
Jn  J  n 


(9a 


Next,  utilize  the  product  rule  for  differentiation,  i.e 
udv  =  d(uv)  -  vdu 


New  the  volume  integrals  of  gradients  may  be  replaced  by  surface 
integrals,  using  the  Guass-Divergence  Theorem. 


£  Wij)  dr  -  JyUi  jGu,.)  dT  +£(\jG).jui  dn 

-  f  ft,-  jGUj )  dT  +  f  (Oj  jGJ  iUj  dn 

Jr  Jn 

dr  +  ^  ■  jn^mQidn=o  (11) 


Reordering  the  terms  and  collecting  under  common  integral  signs: 


j*  nj(Vij)  "  “  ^(dijGUj)  -  M^uK)  dT 

+  +  ^.jG).iuj  +  <V^.A  dn  -  f  u^  dn  =  0 

Jn  n  J 


(12) 


The  indices  can  new  be  rewritten  so  that  the  unknowns,  ,  may  be 
factored  out. 


| _n,<V„)  -  <nA,js+  niflj,iG+  niai.i>  )ui  <3T 

+  f  +  +  <\jx).i]ui  1X1  -  f  *,*■»,  *  =  o  (13) 

Jn  J  n 


njUjj  in  the  first  term  of  Eq.  (13)  is  simply  the  surface  traction  and 
can  be  defined  as  t(  =  nJ-cri j .  For  further  simplification,  the  second 
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term  under  the  surface  integral  will  be  defined  as  t*f  =  n^  -G  +  njUj  (G 
+  r^Uj  jX .  So, 

f  (Uit,  -  t*,u,)dr  +  f  [(^fjG)  -  +  (ajfiG)  -  +  dn 

Jr  Jn 

-  f  dn  =  0  (14) 

Jn 


The  first  domain  integral  of  Eq.  (14)  is  solved  by  setting  the 
assumed  function  terms,  commonly  called  the  adjoint  opertator,  equal  to 
a  vector  of  Dirac  Delta  functions. 


+  (fly®,,  +  (fii.jX),,  =  Sfr.Xj*, 


(15) 


where  6  (x,^)  is  the  Dirac  delta  function  and  §{  is  a  unit  vector. 

The  choice  of  the  Dirac  Delta  function  new  allows  the  domain  integral  to 
be  integrated  exactly.  So,  substituting  Eq.  (15)  into  the  domain 
integral  and  integrating  gives 


dn 


Cu(x) 


(16) 


where  x^  is  the  integration  variable  and  C  is  a  fraction  dependent  upon 
the  location  of  integration.  If  the  sigularity  due  to  the  Dirac  delta 
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function  at  x  =  is  integrated  about  completely,  C  is  equal  to  one. 
For  the  purpose  of  determining  the  boundary  displacements,  the 
integration  will  be  done  on  a  smooth  boundary  surface  and  C  will  be 
equal  to  1/2. 

Making  the  appropriate  substitution,  Eq.  (14)  becomes 


Cu(x)  +  f  (u^.  -  t'jU^dT  -  f  dn  =  0  (17) 

Jr  J  n 

In  order  to  evaluate  the  boundary  integrals,  u,-  must  be  solved  for. 
Ihis  task  amounts  to  solving  Eq  (15) .  In  order  to  solve  Eq.  (15) ,  ul-  is 
written  as  u-  =  ujt  (x,xo)et.  Ihis  second  order  tensor  has  the 
interpretation  that  the  individual  elements  of  uJt  are  the 
displacements  in  the  Jth  direction  at  the  point  due  to  a  unit  point 
force  acting  in  the  Lth  direction,  given  by  et,  applied  at  point  x. 

With  this  understanding,  Eq.  (15)  can  be  written 

(G&jiJ.ii  +  C(G+X  )ujl]>lj  =  & (x,5^j) (18) 

In  2-D,  Danson  has  solved  Eq.  (18)  to  give  (6:211-213) 

=  4JiOH/){(3"4^)lni  5ji  +  rJri>  (19) 


and  therefore 
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t#ji  =  "  47T(l-I/)r  CVka1"2*')^  +  2rjri}  -  (1-2U)  (r^-r^)  ]  (20) 

where  r  is  the  magnitude  of  the  vector  between  the  point  being  solved 
for  (x)  and  the  field  point  (xo) ,  V  is  Poisson's  ratio  and  the  rk's  are 
the  direction  cosines  of  r. 

When  Eqs.  (19)  and  (20)  are  substituted  into  Eq.  (17)  the  surface 
integral  can  be  evaluated  given  a  suitable  displacement  function. 

Except  for  the  mass  term,  which  will  be  handled  latter,  Eq  (17) 
represents  the  boundary  integral  formulation  of  the  elastic-dynamic 
equation. 

2.2  Boundary  Element  Formulation 
Eq.  (17)  can  be  discretized  by  creating  boundary  elements  over  the 
structure.  Each  integral  can  be  written  as  a  summation  of  integrals 
over  each  element: 


Ox 


I  t  "  f  dr  +  f  t*,.^.  dr  ]  +  f  u^-  dn 
k=i  J  rk  J  rk  J  n 


(21) 


For  a  2-D  structure  each  element  would  appear  as  a  line.  Figure  1  shews 
a  simple  2-D  beam  section  paved  with  16  boundary  elements.  By  assuming 
a  shape  function  for  each  element,  the  displacement  at  any  point  on  the 
element  can  be  written  as  a  function  of  the  nodal  point  values.  For 
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instance,  for  a  single  element  as  given  in  Figure  2. 
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Fig.  1.  Beam  Section  with  16  Elements 


a  O - O  b 

Fig.  2.  Single  Element 


the  displacement  at  any  point  x  on  the  element  can  be  written  as 

+  4=£  <22> 

where  L  is  the  length  of  the  element  and  Ug  and  are  the  displacements 
at  the  a  and  b  nodes.  Eq.  (22)  gives  a  linear  relationship  along  the 
element.  Higher  order  elements  can  be  used  (i.e.  quadratic,  etc) 
depending  upon  the  accuracy  requirements.  More  generally,  the 
functional  formulation  over  each  element  may  be  written  as 


u(x)  =  [  Na  Nb  ]{  u*  Uh  } 


(23) 


or 

u(x)  =  [N]{u)  (24) 

where  (u)  is  a  column  vector  of  nodal  displacements,  [N]  is  a  row  of 
shape  functions,  and  u(x)  is  the  value  of  the  displacement  at  point  x  on 
the  element.  Surface  tractions,  t,  can  be  approximated  in  a  similar 
manner. 

With  the  introduction  of  elements  and  shape  functions,  the  boundary 
integral  equation  can  now  be  written  as  a  summation  of  integrals  over 
the  individual  elements  with  displacements  and  tractions  as  functions  of 
the  nodal  value  vectors. 


Cu  = 


1 

I  [" 

k=l 


f  [Ni^t,  dr 

K  r  r 

+  *■  ]  -  I  ® 

Jr  Jn 


mi^u,-  dn 


(25) 


Since  the  nodal  values  under  the  integral  signs  are  constants  they  can 
be  taken  outside  the  integrals.  In  addition,  if  Eq.  (23)  is  written  for 
each  nodal  point,  the  resulting  equations  can  be  recast  in  matrix  form. 


C{u)  =  [G]  (t)  -  (S]{u)  -  (m*) 


(26) 
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where  [S]  and  [G]  are  matrices  of  the  elemental  integrals  associated 
with  each  nodal  point  and  {m* }  are  terms  associated  with  the  mass 
integral  (yet  to  be  determined) .  The  elemental  integrals  can  be 
evaluated  numerically  using  Guassian  quadrature.  The  entire 
discretization  process  is  very  similar  to  that  used  in  finite  element 
methods  and  is  explained  in  depth  by  Gipson  (7:115-120) .  Eq.  (26)  can 
be  further  simplified  by  corrfoining  [S]  and  C: 

[H]  (u)  =  [G] (t)  -  {m*}  (27) 

where  [H]  =  [S]  +  C[I] 

2.3  Treatment  of  the  mass 

Up  until  now  the  mass  integral  has  been  ignored,  but  this 
integration  must  be  completed  in  order  to  solve  the  equation  for  dynamic 
motion.  Ahmad  and  Banerjee  have  developed  a  method  for  handling  the 
mass  term  using  an  approximated  density  function  and  particular 
integrals  (8:682-694).  After  applying  the  method  of  particular 
integrals,  Eq.  (27)  becomes 

[G]  (t)  -  [H]  (u)  =W2([G][Q]  -  [H]  [D] )  [K]  (u)  (28) 

where  w  is  a  frequency  term  resulting  from  the  assumption  of  sinusoidal 
motion  and  the  [Q],  [D],  and  [K]  matrices  are  developed  from  the 
particular  integral  method.  The  [G]  and  [H]  matrices  remain  the  same  as 
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III.  Wave  Propagation 


3.1  Substructure  Analyses 

At  this  point  Eq.  (28)  could  be  applied  to  an  entire  structure.  If 
the  structure  is  repetitive,  such  as  a  truss  structure,  it  is  possible 
to  analyze  only  one  bay  of  the  truss  and  still  calculate  the  dynamical 
behavior  of  the  structure. 

Figure  3  shews  one  bay  of  a  two-dimensional  truss  structure  with 
two  attach  points  on  each  side.  The  tractions  and  displacements 
associated  with  each  connecting  point  sure  also  shewn. 


3  ^3'  tr3 
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Figure  3  Single  Bay  of  a  Repetitive  Truss 


An  equation  relating  the  left  and  right  forces  and  displacements  can 
be  written  in  terms  of  a  transfer  matrix  [T] : 

(vr)  =  [T]{vt}  (29) 


L5 


where 


(vt}  -  (ui3/ut4  ^13  *-14) 


(30) 


<Vr>  » 


(31) 


The  vector  (v)  at  any  junction  on  the  truss  is  called  the  state  vector 
at  that  junction.  Since  each  bay  has  the  identical  structural 
characteristics,  the  transfer  matrix  for  each  bay  will  also  be 
identical.  This  makes  it  possible  to  describe  the  state  at  any  bay 
junction  using  only  the  substructure  transfer  matrix.  For  instance,  the 
forces  and  displacements  at  the  n01  bay  can  be  written  as: 

(Vr}n  =  [T]n{Vl}1  (32) 

Thus,  the  state  vector  at  any  bay  junction  can  be  propagated  along  the 
structure  by  use  of  the  transfer  matrix. 

3.2  Transfer  Matrix  Derivation 

The  transfer  matrix  contains  information  about  the  mass  and  stiffness 
of  the  substructure.  Using  the  BEM  described  in  chapter  2,  the  transfer 
matrix  for  the  substructure  can  be  derived.  Starting  with  Eq.  (28) : 

[G]  (t)-(H]  (U)  =  WZ([G][Q]-[H][D])[K])(U)  (28) 
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The  right  hand  side  of  the  Eq.  (28)  can  be  simplified  into  a  single 
matrix  [m] ,  by  carrying  out  the  appropriate  multiplication  and 
subtraction: 


[m]  =  v^(  [G]  [Q]~[H]  [D] )  [K]  (33) 

Substituting  Eq.  (33)  into  Eq.  (28)  gives 

[G](t}-[H](u)  =  [m]  (U)  (34) 

The  [H]  matrix  can  now  be  added  to  [m]  to  define  a  new  matrix  [M] . 

[M]  =  [m]  +  [H] 

Substituting  Eq.  (35)  into  Eq.  (34)  gives 


[G](t)  =  [M]  (u) 


(36) 


Finally,  multiplying  through  by  [M]"1: 


where 


[R](t)  =  (U) 


(37) 


[R]  =[M]"1  [G] 
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[R]  can  be  thought  of  as  the  dynamical  admittance  matrix  equation  for 
the  boundary  nodes  of  the  substructure. 

Eq.  (37)  can  be  partitioned  to  separate  the  left  connecting,  right 
connecting,  and  outer  surface  nodes  as  follcws. 


Rtl  , 

*lr 

Rio 

Rrr 

Rfo 

•  •  « 


(39) 


Where  t0  and  represent  the  nodal  values  that  do  not  lie  on  a 
junction.  Because  there  is  no  contact  at  the  non- junction  nodes,  the 
traction  on  those  nodes  is  identically  equal  to  zero.  Therefore,  all 
outer  surface  terms  can  be  eliminated  from  Eq.  (39)  to  give 


-  .  - 

-  - 

—  — 

Rll  ;Rlr 

JVl  {Rrr  _ 

5  _ 

_ur_ 

For  simplification,  Eq.  (40  ) 

can  be  written  as 

(40) 
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At  this  point,  seme  care  must  be  taken  to  assure  compatibility 
between  the  sign  convention  of  the  HEM  formulation  and  the  transfer 
matrix  equation.  Figure  4.  shews  the  sign  conventions  vised  for  each 
formulation. 


v  v 


V  V 


Hil'  Si  Hu-'  Sr 

Boundary  Element  Coordinate  Definitions 


Si'  “  Si 

♦ 


Uyr'  Sr 


H(l'  Si  Hu-'  Sr 

Transfer  Matrix  Coordinate  Definitions 


Fig.  4.  Coordinate  System  Definitions 


As  shewn,  the  sign  convention  of  tt  for  the  transfer  matrix  is  opposite 
to  that  of  EEM.  Therefore,  the  left  side  tractions,  in  Eq.  (41)  must  be 
multiplied  by  -1  in  order  to  satisfy  the  transfer  matrix  sign 
convention.  After  multiplying  tt  by  -1  and  rearranging,  Eq.  (41)  can  be 
rewritten  as 


Up 

DB'1  -  C  : 

• 

.......... 

DB_1A 

“l 

Jr_ 

• 

B_1C 

(42) 


or 


(vr)  =  [T](vt) 


(43) 


3*3  Eigenvalues  of  the  Transfer  Matrix 

Using  the  approach  of  Signorelli  (4:23-25),  wave  propagation  in  a 
repetitive  structure  can  be  represented  as: 

(vr)  =  e{vt)  (44) 

Eq  (44)  shews  that  the  state  vector  at  the  right  side  of  the 
substructure  is  the  same  as  the  state  vector  at  the  left  multiplied  by  a 
factor  e.  e  will  generally  be  complex  due  to  the  phase  difference 
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between  the  response  at  each  end  of  the  substructure.  If  Eq  (44)  is 
combined  with  Eq  (43)  the  following  eigenvalue  problem  in  €  is  formed. 

([T]  -  [I]e){Vl>  =  0  (45) 

For  any  given  frequency  in  [T] ,  Eq.  (45)  will  produce  a  set  of 
eigenvalues.  In  addition  to  being  complex,  e  will  occur  in  e  and  1/e 
pairs,  corresponding  to  right  and  left  going  waves.  For  right  going 
waves,  the  magnitude  of  e  will  be  less  than  1.  For  left  going  waves, 
the  magnitude  of  e  will  be  greater  then  1.  Eigenvalue  magnitudes  equal 
to  unity  represent  a  wave  mode  that  will  propagate  undiminished  across 
the  structure.  These  undiminished  wave  modes  are  said  to  be  in  a  pass 
band.  Eigenvalues  with  a  magnitude  other  than  1  represent  non- 
propagating  wave  modes  and  are  said  to  be  in  a  step  band.  Eigenvalue 
behavior  can  be  represented  on  the  e  plane,  as  shewn  in  Figure  5. 


Ime  Ime  Ime 


|e|  <  1 
stop  band 


M  =  1 

pass  band 


I  e  |  >  1 

step  band 


Figure  5. 


Eigenvalues  on  the  e  Plane 
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The  eigenvalues  of  the  transfer  matrix  can  also  be  represented  in 


exponential  form: 


e 


(46) 


Where  L  is  the  length  of  the  bay.  a  is  complex  and  can  be  separated 
into  its  real  and  complex  parts: 


a  =  a  +  bi 


(47) 


so 


€ 


=  e  +  e 


biL 


(48) 


or 

e  =  e*  +  ei<kL  *  2nr>  (49) 

k  is  a  nondimensional  wave  number  and  is  related  to  wavelength  y  by 
k  =  27T/Y-  e31  is  an  attenuation  coefficient  and  describes  the  rate  of 

decay  as  the  wave  passes  through  a  bay.  Negative  values  of  aL  indicate 
a  left  or  negative  going  wave.  The  imaginary  part  of  Eq.  (49) 
describes  the  phase  relationship  between  the  state  vectors  at  the  left 
and  right  side  of  the  bay. 

I 

I 


IV. 


4.1  Wave  Propagation  Of  a  Two-Dimensional  Beam 

A  two-dimensional  beam  was  selected  as  a  baseline  structure  for 
demonstrating  the  implementation  of  boundary  element  analysis  with  wave 
propagation  theory  as  discussed  in  chapters  1,  2,  and  3.  The  selection 
of  a  simple  beam  was  based  on  three  criteria.  First,  results  are 
easily  compared  to  well  understood  continuum  models.  Second,  a  beam 
model  is  simple  to  implement  yet  still  provides  verification  of  the 
1  analysis  procedures.  Third,  a  beam  can  be  thought  of  as  a  repetitive 

structure  if  divided  into  several  short  beam  elements. 

i 

Figure  6  shows  a  long  beam  divided  into  several  sections  or  "bays". 


i 


Figure  6  Sectioned  Beam 


The  properties  of  the  beam  were  selected  as  follows: 

E  =  l.OxlO6  psi 
h  =  5.0  in 
m  =  .10  lb/ in3 
t  =  1.0  in 
V  =  .2 
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where  E  is  the  modulus  of  elasticity,  h  is  height,  m  is  mass  density 
t  is  thickness  and  V  is  poisson's  ratio. 

A  10  inch  length  of  beam  was  selected  to  represent  one  bay.  The 
bay  was  modeled  with  30  linear  boundary  elements  as  shewn  in  Figure  7. 


Figure  7  Boundary  Elements  on  Beam  Section 


A  FORERAN  program  written  by  Brebbia  for  static  2-D  structural  analyses 
was  modified  to  calculate  the  [H]  and  [G]  matrices  of  Eq.  (27)  (9:429- 
438) .  Additional  FORERAN  code  was  written  to  implement  the  particular 
integral  method  discussed  in  section  2.3  (Eq.  (28)).  The  transfer 
matrix  was  then  calculated  using  the  procedure  outlined  in  chapter  3. 

The  6  node  interface  of  the  beam  section  produced  a  24x24  transfer 
matrix  (12  displacements  and  12  tractions) .  Thus,  24  eigenvalues  could 
be  extracted  at  any  given  frequency* 


4.1.1  Eigenvalue  Analysis 


Hie  eigenvalues  of  the  transfer  matrix  were  extracted  for 
frequencies  from  10  to  200Hz  using  an  EISPACK  Fortran  solver  (10:26-27) . 
At  each  frequency,  four  eigenvalues  of  magnitude  1.0  always  appeared. 

As  discussed  in  section  3.3,  eigenvalues  of  magnitude  1.0  indicate  a 
propagating  wave  mode.  Since  complex  eigenvalues  always  appear  in 
complex  conjugate  pairs,  the  four  propagating  eigenvalues  represented 
two  wave  modes. 

Inspection  of  the  eigenvectors  revealed  that  the  two  propagating 
wave  inodes  were  a  bending  mode,  and  a  compression-extension  mode,  as 
expected  for  a  beam.  Figure  8  shews  a  dispersion  curve  of  the  phase  of 
the  eigenvalue  versus  frequency  for  the  bending  mode.  A  similar  plot 
for  the  compression-extension  mode  is  shewn  in  Figure  9.  As  shewn,  the 
phase  angle  increases  with  frequency.  This  is  due  to  the  decreasing 
wavelengths. 

The  remaining  right  going  wave  modes  exhibited  step  band  behavior 
( |  e  |  <  1  )  at  all  frequencies.  A  sampling  of  the  step  band  eigenvalues 
is  presented  in  Table  1.  The  step  band  eigenvalues  represent  localized 
modes  that  quickly  die  out  as  they  travel  down  the  beam.  Eigenvalues 
with  a  magnitude  much  smaller  than  1  represent  localized  behavior  that 
does  not  propagate. 

Classical  beam  theory  provides  a  basis  for  comparing  the  results 
shewn  in  Figures  8  and  9.  For  bending,  simple  2-D  beam  theory  predicts 
the  following  behavior  for  a  free-free  beam  (11:163-166). 
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FREQUENCY  C«z) 


Figure  9 


Compression-Extension  Mode  Dispersion  Curve 


Table  1  Step  Band  Eigenvalues  for  Beam  Bay  at  lOOhz 


lei 

Realfe} 

Imaafe} 

Fhase(dea) 

. 4282E+0 

.4282E+0 

.0 

.0 

.9409E-3 

-.8778E-3 

.  3387E-3 

-.2110E+2 

. 1026E-3 

-.1026E-3 

.0 

.0 

-.6920E-4 

-.2352E-4 

. 6508E-4 

-.7013E+2 

.8253E-6 

.8253E-6 

.0 

.0 

. 3416E-6 

-.3416E-6 

.0 

.0 

W  =  (2 . 57f)2  (EI/mAL4) 1/2  (50) 


where  w  is  the  frequency  for  one  wave  of  a  bending  mode,  I  is  the  moment 
of  inertia,  A  is  the  cross  sectional  area  and  L  is  the  length.  The 
phase  of  the  wave,  Bb,  at  any  point  on  the  beam  can  be  described  by  a 
simple  ratio: 


%  =  2to/L  (51) 

where  x  is  the  distance  of  the  point  frera  the  end  of  the  beam  and  2 it  is 
simply  the  phase  of  one  complete  wave.  Solving  Eq.  (50)  for  1/L 

1/L  =  {(w/(2.57T)2(mA/EI)1/2}1/2  (52) 

and  multiplying  through  by  27 rx 

27TX/L  =27rx((w/(2.5m)2(mA/EI)V2)1/2  (53) 
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gives  the  theoretical  phase  angle: 


Bh  =27rx{(w/(2.57r)2(mVEI)1/2}1/2  (54) 

A  similar  development  for  the  theoretical  compression-extension  phase 
angle  gives. 


Rce=  (xw/ttXiVE)172 


(55) 


If  the  length  of  the  beam  bay  is  used  for  x  in  Eq.  (54)  and  (55) ,  a 
comparison  to  calculated  results  can  be  made.  This  comparison  is  shewn 
in  Figures  10  and  11.  As  shown  in  Figure  10,  there  is  closer  agreement 
with  bending  theory  at  lower  frequencies.  As  the  frequency  increases, 
the  EEM  is  less  able  to  model  the  bending  behavior.  This  disparity  at 
higher  frequencies  is  not  uncommon  for  discreet  formulations  and  is 
similar  to  results  obtained  by  Signorelli  using  FEM  (4:30) .  Seme  of  the 
disparity  may  be  due  to  the  use  of  linear  elements  in  the  formulation. 
Figure  11  shows  nearly  perfect  agreement  with  theory  for  the  compression 
extension  mode.  This  should  be  expected,  because  compression-extension 
deformations  are  more  easily  modeled  with  linear  elements. 

4.1.2  Eigenvector  Mode  Shapes 

For  each  eigenvalue,  a  corresponding  eigenvector  (V)  can  be 
calculated  and  used  to  generate  a  plot  of  the  deformed  structure.  Each 
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FREQUENCY  ^hZ} 

□  TRANSFER  MATRIX 

Figure  10  Theoretical  vs  Analytical  Bending  Mode 
Dispersion  Curve 
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Figure  11  Theoretical  vs  Analytical  Compression-Extension 
Mode  Dispersion  Curve 


eigenvector  contains  values  of  the  displacements  and  forces  on  the  left 
hand  side  of  the  beam  section: 


{vr}  =  {u,,tr}  (56) 

In  turn,  Eq.  (44)  can  be  used  to  calculate  the  state  vector  on  the  right 
hand  side  of  the  section. 


(vr)  =  £{vl}  (57) 

The  right  hand  side  of  one  section  now  becomes  the  left  hand  side  of  the 
adjacent  section,  allowing  the  state  vector  to  be  propagated  along  the 
beam. 


<Vr>n  =  €(Vr>rv1 


(58) 


were  n  is  a  bay  number.  Recovery  of  the  deflections  on  the  outer 
surface  of  any  bay  is  accomplished  by  expanding  Eq.  (39)  to  get 

Hn  =  PWWn  +  [^]{tr)n  (59) 


By  adding  the  real  part  of  the  deflections  to  the  original  node 
locations,  the  mode  shapes  can  plotted. 
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Figure  12  shews  plots  of  the  bending  mode  for  several  frequencies. 
Deflections  are  multiplied  by  an  appropriate  scaling  factor  to 
accentuate  the  mode  shape.  In  addition,  the  plots  of  Figure  12  were 
generated  using  only  the  deflections  of  the  6  junction  nodes  propagated 
to  the  right  and  connected  by  straight  lines  at  the  outer  beam  surface. 
A  closer  look  at  the  bending  behavior  for  a  single  bay,  with  outer 
surface  node  deflections  included,  is  shewn  in  Figure  13. 


40hz  Bending  Mode 


60hz  Bending  Mode 


80hz  Bending  Mode 


lOOhz  Bending  Mode 

Figure  12  Bending  Modes  at  Selected  Frequencies 
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I  Figure  13  Bending  Mode  For  Single 

Bay  at  60hz 

^  The  number  of  bays  necessary  to  complete  one  cycle  of  a  wave  depends 

upon  the  phase  angle  of  the  eigenvalue.  For  instance,  at  80hz  the  phase 
angle  is  48.1  degrees  -  requiring  approximately  7.5  bays  to  complete  one 

|  360  degree  cycle. 

Plots  of  the  compression-extension  mode,  shewn  in  Figure  14,  are  not 
as  dramatic  as  bending,  due  to  the  small  phase  angles  and  relatively 
small  deflections.  As  shewn,  for  200  hz  the  compression-extension  mode 
would  require  approximately  27  bays  to  complete  a  cycle. 

|  i  i  i  i  i  r  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  ~i  i 

Figure  14  Compression-Extension  Mode  at  200hz 

I 

Figure  15  demonstrates  the  behavior  of  a  step  band  mode  at  60hz  with 
e  =  .519  +  O.i.  As  predicted,  the  wave  is  quickly  attenuated.  An 
additional  non  attenuating  wave,  with  e  =  -.92E-3  +  ,21E-3i  is  shown  in 
Figure  16  for  one  bay.  This  localized  wave  does  not  even  have  a 
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noticeable  effect  on  the  right  side  of  the  bay  due  to  the  small  value  of 


It  should  be  noted  that  any  one  wave  mode  does  not  represent  the 
total  behavior  of  a  structure.  The  dynamics  of  any  particular  structure 
would  be  a  time  dependent  linear  combination  of  all  eigenvectors. 


Two  methods  have  been  proposed  for  using  the  transfer  matrix  to 
identify  natural  frequencies.  The  first  method,  proposed  by 
Signerorelli  requires  calculation  of  a  global  transfer  matrix  [T]g 
(4:22).  For  a  structure  with  n  bays,  Eq.  (43)  can  be  applied 


33 


sequentially  to  each  bay  from  left  to  right  to  obtain 


{vr>n  = 


(60) 


so 


[T]g  =  [T]n 


(61) 


and 


<vr)„  =  (Tyv,),  (62) 

Now,  boundary  conditions  can  be  applied.  For  a  cantilever  beam,  the 
boundary  conditions  would  be 


(vr)„  =  «V°>„  <63> 

for  no  tractions  on  the  free  end.  The  secured  end  would  have  zero 
displacements : 


(^l )  i  {0,t^}^ 


(64) 
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So,  after  droping  subscripts  on  the  vectors,  Eq 


(62)  becomes 


{^.,0}  =  [T]g{0,tt} 
[T]g  can  be  partitioned  to  give 


(65) 


(66) 


The  bottom  new  of  Eq.  (66)  gives 

0=[T3){tl)  (67) 

The  only  non-trivial  way  for  Eq.  (67)  to  hold  true  is  for  the 
determinant  of  [T3]  to  equal  zero.  Therefore  a  plot  of  Det(  [T3]  ) 
versus  frequency  should  reveal  the  natural  frequencies.  Unfortunately, 
the  computation  of  [T]n  for  any  significant  number  of  bays  will 
generally  exceed  the  computational  limit  of  the  computer.  An  attempt  at 
computing  det(  [T3]  )  for  just  three  bays  of  the  beam  section  resulted 
in  numerical  overflew. 

Von  Flotcw  proposes  a  second  method  for  determining  natural 
frequencies  (3:516) .  This  method  uses  the  eigenvectors  of  the  transfer 
matrix  along  with  the  boundary  conditions  to  calculate  a  scattering 
matrix  for  each  end  of  the  structure.  By  using  the  scattering  matrices 
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along  with  the  eigenvalues  of  the  transfer  matrix  another  determinant 
problem  is  set  up  by  virtue  of  the  phase  closure  principle.  Again,  this 
method  proves  numerically  unstable  for  any  sizeable  matrices.  The 
application  of  this  method  on  the  beam  section  proved  unsuccessful. 

Although  use  of  BUM  with  wave  propagation  theory  does  not  easily  lend 
itself  to  identifying  natural  frequencies,  there  are  time  domain 
techniques  that  can  be  used  with  wave  theory  to  model  the  real-time 
behavior  of  structures.  Von  Flotcw  has  demonstrated  the  successful  use 
of  a  time  domain  technique  for  the  traveling  wave  analysis  of  a 
structural  network  (3:518). 


V.  Summary 


5.1  Conclusions 

The  purpose  of  this  effort  was  to  demonstrate  the  use  of  boundary 
element  methods  in  performing  a  wave  propagation  analysis  on  periodic 
structures.  The  ccanbination  of  these  two  analysis  techniques  may  be 
useful  in  describing  the  dynamic  behavior  of  repetitive  structures  such 
as  those  planned  as  orbiting  space  platforms.  Past  work  in  the  area  of 
wave  propagation  in  structures  has  been  limited  to  using  continuum  or 
finite  element  structural  models. 

Although  limited  in  its  scope,  this  effort  has  produced  several 
important  results.  First,  a  computer  code  was  produced  that  formulates 
the  transfer  matrix  of  a  2-D  isotropic  element  using  boundary  element 
techniques.  This  code  can  be  used  to  find  the  spatial  eigenvalues  and 
eigenvectors  of  a  periodic  structure.  Second,  results  from  the  analysis 
of  a  baseline  beam  section  were  in  close  agreement  with  theoretical 
results.  The  disparity  of  results  for  bending  at  higher  frequencies  was 
most  likely  due  to  the  use  of  linear  elements  in  the  formulation  of  the 
transfer  matrix.  Disagreement  of  data  at  higher  frequencies  is  not 
unccmmon  in  most  discreet  formulations,  including  FEM.  Third,  although 
BEM  matrices  are  generally  smaller  then  those  of  equivalent  FEM,  their 
dense  nature  can  produce  problems  in  sane  numerical  techniques  such  as 
finding  determinants.  Finally,  the  demonstration  of  this  method  in  2-D 
holds  promise  for  its  use  in  more  ccmplex  structural  descriptions. 
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5.2  Recommendations 


Hie  following  suggestions  are  made  for  future  work  in  this  area. 

(1)  Investigate  the  use  of  EEM  for  describing  the  behavior  of  joints  in 
periodic  structures. 

(2)  Develop  a  code  for  formulating  the  transfer  matrix  of  a  three 
dimensional  structure  using  boundary  element  techniques. 

(3)  Investigate  the  inclusion  of  damping  terms  in  the  boundary  element 
equation. 

(4)  combine  EEM  and  FEM,  where  appropriate,  to  model  large  truss  like 
structures.  For  instance,  FEM  might  best  be  suited  to  rod  and  beam 
elements  while  BEM  could  be  used  for  joints  and  other  solid  body  members 
of  the  structure. 

(5)  Compare  the  computational  efficiency  between  BEM  and  FEM  in 
formulating  the  transfer  matrix  of  different  structures. 

(6)  Use  BEM  in  describing  nonlinear  behavior  such  as  plastic 
deformation. 

(7)  Develop  an  equivalent  transfer  matrix  description  of  nonsymmetric 
2-D  structures  such  as  triangles. 
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Appendix  A:  Ahmad  and  Baneriee's  Particular  Integral  Method 

Ihe  particular  integral  method  uses  the  fact  that  the  solution  to  the 
elastic  dynamic  equations  can  be  written  as  the  sum  of  a  complementary 
function,  uCj,  and  a  particular  integral,  uPj.  Hie  complementary 
solution  satisfies 


GucjMi  +  (C +X  )uci,jl-  =  °  (68) 

while  the  particular  integral  satisfies 

GuPj,  n  +  (G  +X)upj/jj  =  itfoAij  (69) 

Eq.  (69)  still  contains  the  unknown  displacement  field  Uj.  But,  Uj  can 
be  described  using  an  unknown  fictitious  density  function,  f,  and  a 
known  function,  C,  in  a  manner  similar  to  using  shape  functions: 

1 

Uj  -  s  «VV  <7°> 

n=l 

A  simple  function  that  can  be  selected  for  C  is 

Cjk  =  Sik(  R-r)  (71) 


I 
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where  R  is  the  largest  distance  between  two  points  on  the  body  and  r  is 
the  magnitude  of  the  vector  from  the  field  point  to  the  source  point. 
Eq.  (69)  can  now  be  written  as 


1 

GuPj  +  (G  -(X  JuPj  =  mw2^  <S]kfk  (72) 

n=l 

ifj  can  new  be  chosen  as  any  function  that  satisfies  Eq.  (72)  and  can  be 
represented  as 


I  A 

n=l 


(73) 


were  is  found  to  be 


=  mw?  (9-101/ )r  _  11-2V)R 
ik  G  U  90(1-V)  6-8V 


>  V*  ~ 


30(1-1/) 


-jrlcr> 


(74) 


The  surface  traction,  tPjf  related  to  the  displacement  field  of  Eq.  (74) 
is  defined  as 


*j  =  I  <Vic  (75) 

n=l 

Q  can  be  found  using  Eq.  (74)  along  with  the  strain-displacement 
relationships,  and  is  given  by 
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(76) 


Qjk  =  raw2*  ( 


(5L/-11  r 
15(1-1/) 


3-41 ?  rkllj 


+  ( 


(4-5l/)r 
15(1-1 !/) 


(1-2  V)R 
3-4V 


>rjI\ 


+  {( 


(4-5V)r 

15(1-1/) 


(1-2V)R 
3-41 > 


1 

15(1-1/) 


)rini> 


If  Eq.  (27)  is  now  assumed  to  satisfy  the  complementary  solution 

[G](tc)  -[H]{uc>  =  0  (77) 

then 

[G]{t>  -  [H]  (u)  =  (G]{tP>  -  (H]{uP)  (78) 

Substituting  Eq.  (73)  and  (75)  into  (78) ,  gives 

[G]  { t >  -  [H]  (u)  =  ^([G][Q]  -  [H]  [D] )  (f)  (79) 

{f}  can  be  found  from  Eq.  (70) 

{/},•  =  S.j(K]{n)j  (80) 

Where  [K]  is  the  inverse  of  the  symmetric  matrix  formed  by  the  function 
(R-r)  applied  at  each  nodal  point.  Altogether,  Eq.  (79)  becomes 

[G]  (t)  -  [H]{u>  =  W^([G][Q]  -  [H]  [D] )  [K]  (u)  (81) 
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Appendix  B:  Program  for  solution  of  Two  Dimensional 

Transfer  Matrices  by  the  Boundary  Element 
Method 


This  program  reads  an  input  file  that  describes  the  material 
properties  and  geometry  of  the  problem.  Below  is  a  listing  of  the  input 
file  used  for  the  beam  section  of  chapter  4  along  with  an  explanation  of 
the  fields. 


Input  File  Format 


First  line:  Title  of  Problem 

Second  line:  0, (number  of  nodes) , (number  of  elements) , (number  of 

junction  nodes) ,0, (type  of  problem; l=plane  strain  2=plane 
stress) ,0 

Third  line:  (E*thickness) ,  (poisson's 

ratio) ,  (p*thickness) ,  (height) ,  (length) ,  (I) 

Node  position  lines:  (node  number) , (X  coordinate) , (Y 

coordinate) ,  (blank) ,  0 
Note:  The  junction  nodes  must  be  numbered 

consecutively  starting  with  node  1  on  the  right 
side  and  continuing  on  to  the  left  side.  Exterior 
surfaces  are  numbered  last. 

Element  connectivity  lines:  (element  number) , (node  1) , (node  2) 

Note:  The  direction  from  node  1  to  node  2 
must  be  counterclockwise  on  exterior 
surfaces  and  clockwise  on  interior 
surfaces. 


last  line:  0,0 
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RECTANGULAR  BOX  WITH  5X10X1  DIMENSION 

0,30,30,6,0,2,0 

l.Df6, .2,3. 2D-4 ,5,10,10. 041667 

I, 0. ,5. , ,0 
2,0. ,4. , ,0 
3,0. ,3. , ,0 
4,0. ,2. , ,0 
5,0. ,1. , ,0 
6, 0. , 0. , , 0 

7.10. . 5. . .0 

8.10. . 4. . .0 

9.10. . 3. ..0 

10.10. . 2. . .0 

II, 10. ,1.,,0 
12,10. ,0.,,0 

13 .1. . 5. . .0 

14 . 2 . ,  5 . , ,  0 

15.3. . 5. . .0 

16.4. . 5. . .0 

17.5. . 5. . .0 

18.6. . 5. ..0 

19.7. . 5. . .0 

20.8. . 5. . .0 

21.9. . 5. . .0 

22.1. . 0. . .0 

23.2. . 0. . .0 

24.3. . 0. . .0 

25.4. . 0. . .0 

26.5. . 0. . .0 

27.6.. 0...0 

28 .7. . 0. . .0 

29.8. . 0. . .0 

30.9. . 0. . .0 

I, 13,1 

2.14.13 

3.15.14 

4.16.15 

5.17.16 

6.18.17 

7.19.18 

8.20.19 

9.21.20 

10.7.21 

II, 6,22 

12.22.23 

13.23.24 

14.24.25 


15.25.26 

16.26.27 

17.27.28 

13.28.29 

19.29.30 
20,30,12 
21,1,2 

22.2.3 

23.3.4 

24.4.5 

25.5.6 

26.8.7 

27.9.8 

28.10.9 

29.11.10 

30.12.11 

0,0 


Program  Listing 


C  PROGRAM  FOR  CAICOLATION  OF  TOO  DIMENSIONAL 
C  ELASTODYNAMIC  TRANSFER  MATRICES  BY  THE 

C  BOUNDARY  ELEMENT  METHOD 

implicit  real*8(A-H,0-Z) 

COMMON  /RW/  IRE,  IWR 

CCMM3N  /V  D(2,2)  ,XI(6,3)  ,W(6,3)  ,IDUP(50)  ,INC(50,2)  ,C(50) , 
*S(50,3) ,ISYM(100) ,X(100) ,Y(100) ,IFIP(100) ,AK(100,100) ,P(100) , 
*XM(100),AG(100,100) 

*,H4(100, 100) ,  CM  (100, 100) ,TM(100,100) 

INTEGER  CMN 

REAL*8  CM1(100, 100)  ,TM1(100, 100)  ,KM(100,100)  ,WKAREA(2600) 
*,IENIH,  INERT 
C  INH7T 

OPEN(8,FILE='INR7r.IN' , STAIUS= ' OIL ' ) 

OPEN  ( 9 ,  FILE= '  OUTEUT .  OUT ' ,  STAIUS= '  NEW ' ) 

IRE>=8 

IWR=9 

CMN=100 

IDCT=0 

CALL  INIUr(NE,NN,NP, IPL,P0,NN2,NT,NTL,C1,C2, 

*C3 ,  C4 ,  C5 ,  C6 ,  C7 ,  C8 ,  C9 ,  CIO ,  Cll ,  IDSYM,  XSYM,  YSYM,  INFB 
*, R0,E,GM,PMAX,  area,  lenth,  inert) 

C  CCMEUIE  MATRICES  H,G,D,F  AND  K 

CALL  MATRX(NE,NN,NN2,NT,C1,C2,C3,C4,C5,C6,C7,C8, 

*C9 ,  CIO,  Cll ,  PO,  IDSYM,  XSYM,  YSYM,  INFB,  IFA,  NIF,RO,  GM,  EMAX) 

CALL  MAIMUIT  (NN2 , AK,  CM,  CM1) 

CALL  MA3MUIT(NN2,AG,TM,TM1) 

CALL  MATSUB  (NN2 ,  TKL ,  CM1 ,  CM) 
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o  o  o  o  o 


CALL  INVERT  (EM, TM,NN,EMN) 

CALL  MATEXP(NN,TM,KM) 

CALL  MATMULT  (NN2 ,  EM,  KM,  EM) 

CALL  ERESWP  (NN2,NTL,  EM,  AK,  AG,  RO,E,  area,  Denth,  inert 
*,INC,X, Y) 

4  STOP 
END 


1 

C 

2 

C 

3 

61 

60 

4 


SUBROUTINE  INEUT(NE,NN,NP,IEL,P0,NN2,Nr,NTL,Cl,C2, 

*C3  ,C4,C5,C6,C7,C8,C9,  CIO ,  Cll ,  IDSYM,  XSYM,  YSYM,  INFB 
*,  RO,  E ,  GM,KMAX,  AREA,  LENIH ,  INERT) 
implicit  real*8 ( A-H , O-Z ) 

REAL*8  AREA,  LENTH,  INERT 
CCMMDN  /RW/  ERE,  IWR 

COMMON  /A/  D(2,2) ,XI(6,3) ,W(6,3) ,IDUP(50) ,INC(50,2) ,C(50) , 
*S(50,3) ,ISYM(100) ,X(100) ,Y(100) ,IFIP(100) ,AK(100,100) ,P(100) , 
*XM(100) ,AG(100,100) 

*,EM(100, 100) ,EM(100, 100) ,TM(100,100) 

CHARACTER  TmE*70 
WRITE  (IWR,  1) 

FORMAT  (//////////,  24X,  '***  BOUNDARY  ELEMENT  ME 
*THOD  APPLIED  TO  *  *  *'  ,//,24X,  '*  *  *  PLANE  EL 
*ASTOSTATIC  PROBLEMS  *  *  *',///) 

TITLE  OF  PROBLEM 
I?EAD(IRE,2)TTTLE 
FORMAT  (A70) 

WRITE  (IWR,  2)  TITLE 

GENERAL  INFORMATION  ABOUT  THE  PROBLEM 


READ(IRE,  *)  INFB,  NE,NN,NTL,NP,  IPL,  IDSYM 
READ  ( IRE,  *)E,PO,RO,  AREA,  LENIH,  INERT 
FORMAT (II , 14 , 415 , 2F10 . 0) 
IF(INFB.EQ.O)GO  TO  60 
WRITE  (IWR,  61) 

FORMAT (// ,  13X,  ' *  INFINITE  BOUNDARY  *') 
WRITE  (IWR,  4  )NE,NN,NP,  IPL,  IDSYM,  E,P0 


FORMAT (//,15X,  'NO.  ELEMENTS  =' ,  15,//,  15X, 

* 'NO.  NODES  =',I5,//,15X, 'NO.  POINTS  =  ' , 15,//, 15X, 'PROBL. 

*  TYFE=',I5,//,15X,'SYMME.  TYPE  =•  ,15,///, 15X, 'MATERIAL 

♦FHDPERTIES' ,//,  15X,  'E  =',F10.0, 


*//#15X,  'POISSON  =',F15. 5,///, 3 OX,  'COORDINATES  OF  BOUNERY  NODES', 
*//,! 2X,  'NODE' ,  14X,  'X',15X,  'Y',12X,  'DCUBIE',/) 

NN2=NN*2 


NT=NN+NP 

C  NODES  AND  POINTS  COORDINATES 
DO  5  1=1, NN 
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I 

READ(IRE,*)K,X(K)  ,Y(K)  ,IDUP(K)  ,ISYM(K) 

6  P0RMAT(I5,2F10. 0,215) 

IF(IDUP(K)  .EQ.O)GO  TO  5 
J=IDUP(K) 

IDUP(J)=K 
X(K)=X(J) 

Y(K)=Y(J) 

5  CONTINUE 
DO  63  K=1,NN 
IF(IDUP(K)  .NE.O)GO  TO  62 
WRirE(IWR,7)K,X(K)  ,Y(K) 

GO  TO  63 

62  WRITE (IWR,  16) K,X(K)  ,Y(K)  ,IDUP(K) 

16  FORMAT  ( 10X,  15, 5X,  F15 . 4 ,  IX,  F15 . 4 , 7X,  15) 

63  CONTINUE 

7  PORMAT(10X,  15, 5X,  F15 . 4 ,  IX,  F15 .4) 

IF  (NP.EQ.O)GO  TO  9 
WRITE  (IWR,  8) 

8  FORMAT (//,30X,  'COORDINATES  OF  INTERNAL  POINTS'  ,//,llX,  'POINT' , 
*14X,  *X',15X,  'Y',/) 

K=NNfl 

READ(IRE,*)  (J,X(J)  ,Y(J)  ,ISYM(J)  ,JJ=K,NT) 

14  FORMAT  (15 , 2F10 . 0 , 5X,  15) 

WRITE(IWR,7)  (J,X(J)  ,Y(J)  ,J=K,NT) 

C  NODES  AND  POINTS  AT  SYMMETRY  LINES 

9  IF(IDSYM.EQ.O)GO  TO  49 
WRITE  (IWR,  42) 

42  FORMAT (//,30X,  'BOUNERY  NODES  AND  INTERNAL  POINTS  AT  SYMMETRY 
*  LINE'  ,//,12X,  'L.  X',12X, 'L.  Y',/) 

DO  43  K=1,NT 
IF(ISYM(K)  .EQ.O)GO  TO  43 
IZZ=ISYM(K) 

GOTO  (44,45,46) ,IZZ 

44  YSYM=Y(K) 

WRITE  (IWR,  47)  K 

47  FORMAT  (10X,  15) 

GO  TO  43 

45  XSYM=X(K) 

WRITE  (IWR,  48)  K 

48  FORMAT  (26X,  15) 

GO  TO  43 

46  WRITE  (IWR,  50)  K,K 
50  FORMAT (10X, 15, 11X, 15) 

43  CONTINUE 

C  ELEMENT  CONNECTIVITY 

49  WRITE  (IWR,  10) 

10  PQRMAT(//,30X,  'ELEMENT  CONNECTIVITY'  ,//,13X,  'EL'  ,13X,  'N.  1', 
*12X, 'N.  2' ,14X, 'L' ,/) 

DO  11  1=1, NE 

R£AD(IRE,*)K, INC(K,  1)  ,INC(K,2) 
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12  P0RMAT(3I5) 

II=INC(K,1) 

IF=INC(K,  2) 

11  C(K)=SQRT(  (X(IF)  -X(II) )  **2+(Y(IF)  -Y(II) )  **2) 

WRITE(IWR,13)  (1,010(1,1)  ,INC(I,2)  ,C(I)  ,I=1,NE) 

13  FORMAT (10X, 15 , 11X, 15 , 11X, 15, 5X, F15 . 4) 

C  CONSTANTS 

GH=E/ (2.*(1.+P0) ) 

C11=F0 

IF(IPL-1)40,40,41 

40  PO=PO/(1.+PO) 
dl=0. 

41  C2=3.-4.*PO 

C3=l ./ ( ( 1 . -PO) *12 . 56637062 ) 

C4=l.-2.*PO 

C6=2.*C3*GM 

C7=l.-4*FO 

C1=C3/(2.*GM) 

C5=Cl/2. 

C8=2.*GfV(l.-P0) 

C9=P0/(1.-P0) 

C10= ( 2 . -PO)  /  ( 1 . -PO) 

C  BOUNERY  VAIUES  PRESCRIBED 
DO  19  1=1,  NN2 
P(I)=0 

19  IFIP(I)=0 

READ  ( IRE,  *)NFIP,NDFIP 

20  F0RMAT(2I5) 

WRITE  (IWR,  21)  NFIP,NDFIP 

21  F0RMAT(//,15X,'N0.  DISPL.  FRESC.  =• ,I5,//,15X, 'NO.  TRACT.  ERESC. 
*  =',I5,///,15X, 'DISPIACEMENIS 

*',//,  12X,  'NODE'  ,14X,  'U',15X,  'V',/) 

IF(NFIP.EQ.O. )G0  TO  22 
DO  23  I=1,NFIP 

READ(IRE,*)K,P(2*K-1)  ,P(2*K)  ,IFIP(2*K-1)  ,IFIP(2*K) 

24  FORMAT (15, 2F10. 0,215) 

INI>=IFIP(2*K-1)+2*IFIP(2*K) 

GO  TO  (25,26,27) ,IND 

25  WRITE  (IWR,  28)  K,  P  (2*K-1) 

28  P0RMAT(10X,I5,5X,F15.4) 

GO  TO  23 

26  WRITE  (IWR,  29)  K,P(2*K) 

29  F0EMAT(10X,I5,21X,F15.4) 

GO  TO  23 

27  WRITE(IWR,30)K,P(2*K-1)  ,P(2*K) 

30  PQFMAT(10X,I5,5X,F15.4, 1X,F15.4) 

23  CONTINUE 

22  IF(NDFIP.EQ.O)GO  TO  31 
WRITE  (IWR,  34) 

34  FORMAT (//  ,15X,  'TRACTIONS' ,//,  12X,  'N0DE',13X,  'PX'14X,  'P Y',/) 
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DO  32  1=1 , NDFIP 
READ(IRE,*)K,P(2*K-1) ,P(2*K) 

33  P0FMAT(I5, 2F10.0) 

32  WRITE(IWR,  30)K,P(2*K-1)  ,P(2*K) 

C  OCMPUTE  LARGEST  DISTANCE  BETWEEN  POINTS 

31  EMAX  =0. 

DO  90  1=1, NN 
DO  90  J=1,NN 
XX=X(I)-X(J) 

YY=Y(I)-Y(J) 

EMA=SQRT  (XX**2+YY**2 ) 

IF(EMA.GE.PMAX)  FMAX=PMA 
90  CONTINUE 
C  INTEGRATION  POINTS 

XI  ( 1 , 3) =-0 . 932469514203152 
XI (2 , 3 )  =-0 . 661209386466265 
XI (3 , 3) =-0. 238619186083197 
XI(4,3)=-XI(3,3) 

XI  (5, 3)— XI  (2, 3) 

XI(6,3)=-XI(1,3) 

W(l, 3) =0. 171324492379170 
W  (2, 3) =0.360761573048139 
W(3 , 3) =0 . 467913934572691 
W(4,3)=W(3,3) 

W(5,3)=W(2,3) 

W(6,3)=W(1,3) 

XI  ( 1 , 2 )  =-0 . 86113 63 11594 053 
XI  (2 , 2)  =-0. 339981043584856 
XI(3,2)=-XI(2,2) 

XI(4,2)=-XI(l/2) 

W  ( 1 , 2 ) =0 . 347854845137454 
W(2 , 2 ) =0 . 652145154862546 
W(3,2)=W(2,2) 

W(4,2)=W(1,2) 

XI (1 , 1) =-0 . 577350269189626 
XI(2,1)=-XI(1,1) 

W(l,l)=l. 

W(2, 1)=1. 

REIURN 

END 

C 

C 

C 

C 

C 

C 

C 

SUBROUTINE  MATRX(NE,NN,NN2,NT,C1,C2,C3/C4,C5,C6,C7,CS, 
*C9 ,  CIO ,  Cll ,  PO,  IDSYM,  XSYM,  YSYM,  INFB,  IFA,  NIF,  RO,  GM,  EMAX) 
inplicit  real*8(A-H,0-Z) 
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OCMMON  /A/  D(2,2)  ,XI(6,3)  ,W(6,3)  ,IDUP(50)  ,INC(50,2)  ,C(50) , 
*S(50,3) ,ISYM(100) ,X(100) ,Y(100) ,IFIP(100) ,AK(100,100) ,P(100) , 
*XM(100)  ,AG(100, 100) 

*,1*1(100,100) , CM (100, 100)  ,1m  (100, 100) 

REAL*8  NORM1 (100)  ,NQRM2 (100)  ,NOPM(2) 

OCMMON  /A4/  H(3,4)  ,G(3,4)  ,HL(3,4)  ,GL(3,4)  ,NORMl,NOFM2 
*,NQEM,FM(2) 

C  KRONECKER  DELIA 

D(l,l)=l. 

D(2,2)=l. 

D(l,2)=0. 

D(2,l)=0. 

C  CLEAR  ARRAYS 

DO  1  1=1,  NN2 
XM(I)=0. 

DO  1  J=1,NN2 
AK(I,J)=0. 

1  AG(I,J)=0. 

C  OCMTOTE  PARAMETERS  FOR  SYMMETRY  LOOP 
IFAr=l 
NIF=1 

IF  (IDSYM.  EQ.  1)  IFA=2 
IF(IDSYM.NE.2)GO  TO  60 
IFA=3 
NIF=2 

60  IF(IDSYM.EQ.3)IFA=4 
C  TEST  FOR  INFINITE  BOUNERY 
C  IF(INFB.EQ.O)GO  TO  90 
C  DO  91  1=1, NN2 
c  IF(IFIP(I) .NE.O)GO  TO  92 
C  A(I,I)=1 
c  GO  TO  91 
C  92  XM(I)=-P(I) 

C  91  CONTINUE 
C  SYMMETRY  LOOP 
90  DO  2  ISY=1 , IFA , NIF 

C  COMEUTE  CHANGE  SIQi  CCX7TRDLLING  PARAMETERS 
GOTO  (70,71,71,73) ,ISY 
71  IIS=4-ISY 
IFS=IIS 
GO  TO  70 
73  IIS=1 
IFS=2 

C  IOOP  OVER  BOUNERY  NODES 
70  DO  2  1=1, NN 
XS=X(I) 

YS=Y(I) 

IF ( ISY . EQ . 2 . OR.  ISY . EQ . 4 ) YS=2 . * YSYM-YS 
IF(ISY.GE.  3)  XS=2 .  *XSYM-XS 
C  GENERATE  MATRIX  H  AND  G 
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C  NOTE:  OOIDMNS  OF  A  THAT  CORRESPOND  TO  G  ARE  MUITIPIED 
C  BY  C8  TO  AVOID  ROUNDOFF  ERRORS  IN  THE  SOLUTION 
DO  10  J=1,NE 
II=INC(J, 1) 

IF=INC  (J,2) 

IOOD^l 

IF(ISY.NE.1.AND.ISYM(I)  .NE.  (ISY-1)  )G0  TO  6 
IF ( I . EQ . II . OR.  I . EQ . IDUP ( II ) ) I00D=2 
IF(I.EQ.IF.OR.I.EQ.IDUP(IF) )  I00D=3 
6  CALL  FUNC(IGOD, J,C1,C2,C3,C4,C5,C6,C7,P0,II,IF,XS, YS,ISY,IIS, 

*IFS) 

DO  10  K=l,2 
JJ=2* (I — 1) +K 
M=0 

DO  10  NX=1,2 
DO  10  NV=L,2 
M=Mfl 

10=2 *INC ( J,NX) +NV-2 
C  IF(IFIP(IC)  .NE.OJGO  TO  67 
AK(  JJ,  IC)  =AK(  JJ,  IC)  +H(K,M) 

C  XM(JJ)=XM(JJ)4G(K,M)*P(IC) 
c  GO  TO  68 

67  AG(JJ,IC)=AG(JJ,IC)4G(K,M) 
c  XM(JJ)  =XM(JJ)  -H(K,M)  *P(IC) 

C  COMPUTE  REMAINING  OOEFFICIENIS  BY  APPLYING  RIGID  BODY  TRANSIATIONS 

68  GO  TO  (61, 62, 63,64), ISY 

62  IF (NV-2) 61,64,61 

63  IF  (NV-1) 61,64, 61 

64  H(K,M)=~H(K,M) 

C  61  IF  ( IFIP  ( JJ+NV-K)  .  NE .  0 )  GO  TO  69 

61  AK  ( JJ ,  JJ+NV-K)  =AK  ( JJ ,  JJ+NV-K)  -H(K,M) 

GO  TO  10 

C  69  XM(JJ)=XM(JJ)+H(K,M)*P(JJ+NV-K) 

10  CONTINUE 
2  CONTINUE 

C  OCMEUTE  MASS  MATRIX  USING  PARTICULAR  INTEGRALS 
CM1=(9-10*PO)/(90-90*PO) 

CM2= ( 1-2 *P0) / ( 6-8 *P0) 

CM3=1/(30-30*PO) 

CM4=  (  5*P0-1)  /  ( 15*  ( 1-FO) ) 

CM5=2*PO/(3-6*PO) 

CM6= (4-5*P0) / ( 15* ( 1-PO) ) 

DO  300  1=1, NN 
DO  305  J=1,NN 
FM(1)=X(I)-X(J) 

RM(2)=Y(I)-Y(J) 

N0PM(1)=NOTM1(I) 

NCRM(2)=N0FM2(I) 

R2=SQKT(RM(1)  **2+FM(2)  **2) 

EM(I,J)  =  (EMAX-R2) 
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DO  310  K=l, 2 
DO  315  L=l,  2 
IM=2*I+K-2 
JM=2*J+L-2 

c  COMPUTE  D  MATRIX 

EM(IM,JM)=RQ*(  (CM1*R2-CM2*EMAX)  *D(K,L)  *R2**2- 
*  (CM3*PM(K)  *EM(L)  *R2) ) /GM 
IF (R2. IE.. 001) GO  TO  330 
C  COMPUTE  T  MATRIX 

TM(IM,JM)  =RO*  ( (CM4*R2-CM5*EMAX)  *RM(L)  *NORM(K)  +  (CM6*R2-2*CM2*EMAX) 
**FM(K)  *NOIM(L)  +  (  (O16*R2-2*CM2*10. )  *D(K,L)-(2*CM3*RM(K)  *RM(L)  )/R2) 
**(RM(K)  *NOKM(K)+RM(L)  *N0EM(L) ) ) 

GO  TO  315 

330  TM(IMfJM)  =FO*  (  (CM4*R2-CM5*RMAX)  *RM(L)  *N0RM(K)  +  (CM6*R2-2*CM2*RMAX) 
**RM(K)  *NOEM(L) ) 

315  CONTINUE 
310  CONTINUE 
305  CONTINUE 
300  CONTINUE 
RETURN 
END 
C 
C 
C 
C 
C 
C 
C 
C 

SUBROUTINE  FUNC(IOOD,JA,C1,C2,C3,C4,C5,C6,C7,PO,II,IF,XS,YS,ISY, 
*IIS,IFS) 

C  INTEGRAIS  OVER  BOUNERY  ELEMENTS 
implicit  real*8(A-H,0-Z) 

COMMON  /A/  D(2,2) ,XI(6,3) ,W(6,3) ,IDUP(50) ,INC(50,2) ,C(50) , 
*S(50,3) ,ISYM(100) ,X(100) ,Y(100) ,IFIP(100)  ,AK(100,100) ,P(100) , 
*XM(100) ,  AG  (100 ,100) 

*,FM(100,100)  ,EM(100, 100)  ,TM(100,100) 

REAL*8  NOEMl(lOO)  ,NORM2(100)  ,NORM(2) 

COMMON  /A4/  H(3,4)  ,G(3,4)  ,HL(3,4)  ,GL(3,4)  ,N0FM1,N0RM2 
*,N0RM,FM(2) 

DIMENSION  DXY(2)  ,EN(2)  ,B(2)  ,CR(2)  ,UL(2,2)  ,PL(2,2)  ,ULL(2,2,2) , 
*PLL(2,2,2) 

DO  5  KK=1, 3 
DO  5  D=l,4 
GL(KK,L)=0. 

HL(KK,L)=0. 

G(KK,L)=0. 

5  H(KK,L)=0. 

DXY(1)=X(IF)-X(II) 

DXY(2) =Y (IF) -Y (II) 
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BN(1)=DXY(2)/C(JA) 

BN(2)=-DXY  (1)/C(JA) 

N0RM1(IF)=BN(1) 

N0FM2 (IF)=EN(2) 

N0EM1(II)=BN(1) 

N0EM2(II)=BN(2) 

GO  TO  (1,2, 2,1) ,ICOD 
C  SELECT  NO.  OF  INTEGRATION  POINTS 

1  NhS’ 5*SQRT(  <2 ‘ *XS"X(II)  ~X(IF) )  **2+  (2 ‘ *YS-Y(II) -Y(IF) )  **2) /C(JA) 

IF  (SEL.  LE.  1 . 5)  NPI=6 
IF  (SEL.GT.  5 . 5)  NPI=2 
INP=NPI/2 

C  COMPOTE  MATRICES  NUMERICALLY 

IX)  50  KK=1,NPI 

XMXI=0 . 5*  ( 1 .  +XI  (KK,  INP) )  *DXY  ( 1)  +X  ( II)  -XS 
YMYI=0 . 5*  (l.+XI  (KK, INP) )  *DXY (2) +Y (II)  -YS 
R=SQRT(XMXI**2+YMYI**2) 

B  ( 1)  =-0 . 25*  (XI  (KK,  INP)  -1 . )  *C  ( JA) 

B (2 )  =0 . 25*  (XI  (KK,  INP)  +1 . )  *C ( JA) 

ER(1)=XMXI/R 

ER(2)=YMYI/R 

DR0N=DR(1)  *BN(1)+ER(2)  *BN(2) 

C  OCMEOTE  MATRICES  H  AND  G 

DO  6  1=1,2 
DO  6  J=l,2 

UL(I,  J)=-Cl*(C2*DIOG(R)  *D(I,  J)  -ER(I)  *ER(J)  ) 

PL(I,  J)  =-C3*  (  (C4*D(I,  J)+2 .  *DR(I)  *ER(J) )  *EREN+C4*  (ER(J)  *BN(I)  -ER 

6  continue 

DO  7  IA=1, 2 
10=0 

DO  7  LD=1,2 
DO  7  JJ=1,2 
IOIC+1 

G  (LA,  IC)  =G  (LA,  IC)  -KJL(LA,  JJ)  *B(LL)  *W(KK,  INP) 
H(IA,IC)=«(1A,IC)+PL(IA,JJ)  *B(LL)*W(KK,INP) 

7  CONTINUE 

C  OCMEOTE  MATRICES  HL  AND  GL  (INTERNAL  STRESSES) 

10  DO  11  1=1,2 
DO  11  J=I,2 
DO  11  K=l,2 

ULL(I,J,K)=C3*(C4*(CR(J)  *D(K,I)+ER(I)*D(K,J)-CR(K)*D(I,J)+2.*CR(I) 
*)*ER(J)*DR(K))/R  1  ; 

Bl=2.  *EREN*  (C4*ER(K)  *D(I,  J)+PO*  (ER(J)  *D(I,K)+CR(I)  *D(J,K)  -4.  *ER(I) 
*)  *ER(J)  *131  (K) ) 

B2=2 .  *P0*  (EN ( I)  *CR ( J)  *CR (K)  +EN ( J)  *CR ( I)  *ER (K)  ) 

B3=C4*(2.*BN(K)  *CR(I)  *ER(J)+EN(J)  *D(I,K)+EN(I)  *D(J,K) ) 

11  PLL(I,J,K)=C6*(B1+B2+B3-C7*BN(K)*D(I,J))/R**2 
IL=0 
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DO  12  1=1,2 
DO  12  J=I,2 
IL=ILfl 
IC=0 

DO  12  IAA=1,2 
DO  12  JAA=1,2 
IC=IC+1 

GL(IL,IC)  =<3L ( IL,  IC)  +B ( IAA)  *ULL ( I ,  J ,  JAA)  *W (KK,  INP) 

12  HL(IL,  IC)=fttj(IL,  IC)  +B(IAA)  *PLL(I,  J,  JAA)*W(KK,  INP) 

50  CONTINUE 
GO  TO  18 

C  OCMFUTE  MATRICES  H  AND  G  ANALYTICALLY  (BOUNERY  CONSTRAINT  EQ.) 
2  AD=C5*C2*C(JA) 

AAfAL*  (  0 . 5-DIOG  (C  ( JA) ) ) 

DO  15  1=1,2 
DO  15  J=l,4 
IT=(J/2)*2+2-J 

G(I,  J)=C5*DXY(I)  *DXY(IT)/C(JA) 

IF(IT.EQ.  I)  G  (I ,  J)  =G  (I ,  J) +AA 
15  CONTINUE 
IAA?=-2 

IF(I00D.EQ.3)  IAA=0 
G  (1, 3+IAA)  =G(1, 3+IAA) +AL 
G  (2 , 4+IAA)  =G  (2 , 4+IAA)  +AL 
H  ( 1 , 2 -IAA)  =C3*C4*  ( 1 .  +IAA) 

H  ( 2 , 1-IAA)  =-H  ( 1 , 2-IAA) 

C  SYMMETRY  TEST 

18  IF(ISY.EQ.l)GO  TO  8 
DO  24  I=IIS , IFS 
DO  24  J=l,4 
H(I,J)=-H(I,J) 

24  G(I,  J)=-G(I,  J) 

IF(IC0D.NE.4.OR.ISY.EQ.4)GO  TO  8 
DO  25  J=l,4 

HL(2,  J)=-HL(2,  J) 

25  GL(2,  J)=-GL(2,  J) 

8  RETURN 

END 


ADD  TWO  MATRICES 
SUBROUTINE  MATADD(N,A1,A2,A3) 
implicit  real*8(A-H,0-Z) 

REAL*8  Al(100,100) ,A2(100,100)  ,A3(100,100) 
DO  10  1=1, N 
DO  10  J=1,N 

A3  (I ,  J)  =A1  (I ,  J) +A2  (I ,  J) 
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10  CONTINUE 

REIURN 
END 
C 
C 

c 

C  MULTIPLY  TWO  MATRICES 

SUBROUTINE  MATMUIT(N,M1,M2,M3) 
implicit  real*8(A-H,0-Z) 

REAL*8  Ml (100, 100) ,M2(100,100) ,M3(100,100) 

DO  5  1=1, N 
DO  5  J=1,N 
M3 (I, J)=0. 

5  CONTINUE 
DO  10  1=1, N 
DO  10  J=1,N 
DO  10  K=1,N 

M3  (J,  I)  =M3  ( J,  I)  4M1  (J,K)  *M2  (K,  I) 

10  CONTINUE 
REIURN 
END 

C  MULTIPLY  TWO  COMPLEX  MATRICES 

SUBROUTINE  MATMUITI(N,M1,M2,M3) 
implicit  real*8(A-H,0-Z) 

OCMPLEX*16  Ml (100, 100)  ,M2  (100, 100) , M3  (100, 100) 
DO  5  1=1,  N 
DO  5  J=1,N 
M3  (I, J)=(0.D0,0.D0) 

5  CONTINUE 

DO  10  1=1, N 
DO  10  J=1,N 
DO  10  K=1,N 

M3  ( J,  I)  =M3  (J,  I)  4M1  ( J,  K)  *M2  (K,  I) 

10  CONTINUE 

REIURN 
END 


Formulates  the  [K]  matrix  of  the  particular  integral  method 
SUBROUTINE  MATEXP(N,EX1,EX2) 
implicit  real*8(A-H,0-Z) 

REAL*8  EX1(100,100) ,EX2 (100,100) ,D(2,2) 

D(l, 1)=1. 

D(2,2)=l. 

D(l,2)=0. 

D(2,l)=0. 

DO  10  1=1, N 
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DO  10  J=1,N 
DO  10  K=l,2 
DO  10  L=l,2 
II=2*I+K-2 
JJ=2*J+L-2 

EX2  (II, JJ)=EX1(I,  J)  *D(K,L) 

io  continue 

REIURN 

END 

C 

C 

C 

C  SUBTRACT  TWO  MATRICES 

SUBROUTINE  MATSUB(N,SUB1,SUB2,SUB3) 
implicit  real*8  (A-H,0-Z) 

DIMENSION  SUB1(100, 100) ,SUB2 (100, 100) ,SUB3 (100, 100) 

DO  10  1=1, N 
DO  10  J=1,N 

SUB3  (I ,  J)  =SUB1  (I,  J) -SUB2  (I ,  J) 

10  CONTINUE 

RETURN 
END 
C 
C 

c 

c 

c 

SUBROUTINE  FRESWP (NN2,NTL, Ml, KL,G,RO,E, AREA, LENIH, INERT 
*,INC,XO,YO) 

c  This  subprogram  multiplies  the  mass,  [m]  of  Eq  33,  matrix  by  a 

c  given  frequency,  and 

c  then  forms  the  transfer 

c  matrix.  Once  the  transfer  matrix  is  formed  it  is  input  into  the  RGO 
c  routines  for  eigenvalue  and  eigenvector  extraction, 

c 

c  On  input: 

c 

c  KL  is  the  matrix  formulated  in  the  boundary  element  calculations  that 

c  multiplies  the  displacement  terms  (  this  is  the  H  matrix  ) 

c 

c  Ml  is  the  mass  matrix  calculated  by  the  particular  integral 
c  technique, 

c 

c  G  is  the  matrix  formulated  in  the  boundary  element  calculations  that 
c  multiplies  the  displacement  terms, 

c 

c  RO  and  E  are  material  properties.  NN2  is  the  size  of  the  above 
c  matrices. 

c  MEL  is  the  number  of  nodes  on  the  interface  of  the  substructure, 

c  XO, YO,  INC  are  geometry  descriptions  used  for  post  processing 
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IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  Ml(100, 100) ,10.(100,100) ,G(100,100) ,GG(100,100) 

*,MM(100,100) , IM(100, 100)  ,WKAREA(2600)  ,M2(100,100)  ,WR(100) 

*,WI(100)  ,ZR(100,100)  ,FV1(100) , inert, lentil,  IMAGX  (100) 

*,ZI(100,100) ,  IMAGY(IOO) ,REALX(100) ,KEAL¥(100) ,G1(100,100) 
*,01(100,100)  ,EIMAG(100)  ,PART5(100,100)  ,PART6(100, 100) 

*,FLR(50) ,FLI  (50) ,FPR(50) ,FRI(50) ,XO(100) ,Y0(100) 

INTEGER  EMN,IV1(100)  ,N,NTL4, INDEX (100)  ,INC(50,2) 
OPEN(10,FIIE='VECT0R.0OT' ,STAIUS='NEW' ) 

EMN=100 

IDGT=0 

NTL4=4*NTL 

C  this  loop  changes  the  frequency  term  if  desired 
DO  100  N=l, 1 
WRITE (9, 110) 

110  FORMAT  (////) 

HERTZ=  100. 

fre=hertz* (2*3 . 141592654) 
c 

c  thphi  is  the  theoretical  value  of  the  phase  angle  for  ccampresional- 
c  exterrtional  motion  of  a  beam, 

c 

thphi=dsqrt  ( ro/e)  *hertz*lenth*360 
c 

c  thphi2  is  the  theoretical  value  of  the  phase  angle  for  bending  of  a 
c  beam 

c 

thphi2=dsqrt  ( fre*dsqrt  (ro*area/  (e*inert) ) ) 
thphi2=thphi2*lenth*540/  (2.5*3. 141592654) 

DO  20  K=1,NN2 
DO  20  L=1,NN2 
M2  (K,  L)  =  (FEE**2 )  *M1  (K,  L) 

20  OOMTINUE 

CALL  MATADD(NN2,KL,M2,MM) 

CALL  INVERT  (MM,  M2,  NN2,CMN) 

CALL  MA3MJIT(NN2,M2,G,GG) 

C - -C — 


C  MATPART  takes  the  dynamical  admittance  matrix  G  and  manipulates  it 
c  into  the  substructure  transfer  matrix  GG.  Ihis  is  done  by 
c  paritioning 

c  the  matrix,  eliminating  the  interior  nodes  and  rearanging  terms  in 
c  order  to  seperate  the  left  and  right  nodes, 

c 

CALL  MATPART  (GG,NIL,  PART5 ,  PART6 ,  NN2) 
matz=l 

c******************************************** 

c  RGO  calls  the  EISPACK  routine  for  extracting  the  eigenvalues  and 
c  eigenvectors  of  the  transfer  matrix.  The  EISPACK  routines  are  not 


56 


c  included  in  this  listing. 

c******************************************** 

c 

CALL  RG0(EMN,NTL4,GG,WR,WI,MATZ,ZR,ZI,IV1,FV1,IERR) 

c  The  rest  of  the  FRESWP  is  merely  manipulation  of  the  Eigenvalues  and 
c  eigenvectors  for  output, 

do  30  I=l,ntl4 

eimag(I)=dsqrt(wr(I)  **24wi(I)**2) 

C  IF (TEST. NE. 0.0)  1HEN 

C  TEST=0. 

C  GOTO  30 

C  ENDIF 

C 

c  eimag  is  the  magnitude  of  the  eigenvector 
c 

c  phee  is  the  phase  angle 

c 

phee=datan(wi(I)/wr(I) ) *(360.0/6.2831850) 

TEST  =  WI(I) 

TESTMAG=DABS  (EIMAG(I)  -1) 

C  IF  (EIMAG(I)  .gt.l.lODO)  GOTO  30 

C  IF(EIMAG(I)  .LT.l.D-4)  GOTO  30 

WRITE (9, 50) 

50  FORMAT  (/,2X' FREQUENCY  (hertz)  •  ,5X,  'OCMFR  PHASE  (deg)  '  ,5x, 

*'BEND  PHASE  (deg)  ') 

WRITE  ( 9 , 55 )  hertz ,  thphi ,  thphi2 
55  FQRMAT(E15.5/5X,E15.5/5X,  E15.5,/) 

WRITE (9, 60) 

60  FORMAT  (6X' REAL'  ,  13X,  '  IMAGINARY 1 ,  7X,  'EIMAG',  10X,  'PHASE') 

TEST=WI(I) 

WRirE(9,65)WR(I)  ,WI(I)  ,EIMAG(I)  ,FHEE 
65  PCRMAT(E15.5,2X,E15.5,2X,E15.5,2X,E15.5,/) 

WRITE (9, 70) 

70  FORMAT (IX 'NODE  #' ,18X, 'XDISP' ,31X, 'YDISP') 

WRITE (9, 75) 

75  FORMAT (16X,  'REAL',10X,  'IMAG',16x,  'REAL',16X,  'IMAG') 

DO  15  J=l,ntl 

IX=<J*2-1 

IY=J*2 

IMAGX(j)=ZI(IX,I) 

IMAGY(j)=ZI(IY,I) 

KEALX(J)=ZR(IX,I) 

FEALY(J)=ZR(IY,I) 

10  WRITE(9,40)J,ZR(IX,I)  ,IMAGX(J)  ,ZR(IY,I)  ,IMAGY(J) 

40  format (15 , 4X, E15 . 5 , 2X, E15 . 5 , 4x, E15 . 5 ,  lx, E15 . 5) 

15  continue 

DO  150  J=1,2*NTL 
R=2*NTLKJ 
FLR(J)=ZR(K,  I) 
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FLI(J)=ZI(K,I) 

FRR(J)=ZR(K,  I)  *WR(I)  -ZI  (K,  I)  *WI  (I) 

FRI  ( J)  =ZR(K,I)  *WI  (I)  +ZI  (K,  I)  *WR(  I) 

150  CONTINUE 

DO  160  K=l,NN2/2-2*NTL 

KL=NTL*2+K 

IX=2*K-1 

IY=2*K 

IMAGX(KI)=O.DO 
IMM3Y  (KI)  =0 .  DO 
REALX(KE)  =0.  DO 
REALY  (KI)  =0.  DO 
DO  170  J=1,2*NTL 

IMAGX(KI)  =-FLI  (J)  *PART5  (IX,  J)  +FRI  (J)  *PART6  (IX,  J)  +IMAGX(KI) 
IMAGY(KI)=-FLI(J)  *PART5  (IY, J)+FRI  (J)  *PART6  (IY,  J)+IMAGY(KI) 
REALX(KI)=-FLR(J)  *PART5(IX,J)+FRR(J)  *PART6(IX,  J)+REALX(KE) 
REALY (KT)  =-FIR(J)  *PAPT5  ( IY ,  J)  +FRR ( J)  *PAKT6(IY,  J)+REALY(KT) 
170  CONTINUE 
160  OCmTNUE 

c  Vector  is  a  laser  supplied  routine  for  doing  structure  plots 
C  CALL  VECTOR  (NTL,REAIX,REAL¥,IMAGX,IMAGY,WR(I)  ,WI(I)  ,INC 

C  *,XO, YO,NN2/2) 

30  continue 

100  CONTINUE 
return 
end 
C 
C 
C 
C 
C 

SUBROUTINE  MATPART  (GG ,  NTL,  PART5 ,  PART6 ,  NN2 ) 

IMPLICIT  REAL*8  ( A-H ,  O-Z ) 

DIMENSION  PARTI (100, 100) ,PAPT3  (100,100) ,PART4(100,100) , 
*A(100, 100) ,B(100, 100) ,C(100,100) ,D(100,100) , 

*BINV(100,100) ,GG(100, 100) ,PART5(100,100) ,PART6 (100, 100) 
INTEGER  EMN 
EMN=100 

DO  10  I=1,2*NTL 

DO  10  J=1,2*NIL 

A(I,J)=GG(I,J) 

k=2*ntl+j 

B(I,J)=GG(I,k) 

l=2*ntl+i 

C(I,J)=GG(1,J) 

D(I,J)=<3G(l,k) 

10  CONTINUE 

DO  50  1=1 ,  NN2  -4  *NTL 

R=4*NTI>I 

DO  50  J=1,2*NTL 
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JI=2*NTLhJ 

PART5(I,J)=GG(K,J) 

PART6(I,J)=CG(K,JI) 

CONTINUE 

ntl2=2*ntl 

CALL  INVERT (B,BINV,Nri2,EMN) 

CALL  MATMDIT(NTI2,BINV,A,PART4) 

CALL  MATMUIT  (NT12 ,  D,  BINV,  PARTI) 

CALL  MA3MJLT  (NTL2 ,  D,  PART4 ,  PART3 ) 

CALL  MATSUB  (NTL2 ,  PART3 ,  C,  PART3 ) 

DO  30  I=1,2*NTL 

DO  30  J=1,2*NTL 

GG(I,J)=PART1(I,J) 

k=2*ntl+j 

l=2*ntl+i 

GG(I,k)=PART3  (I,  J) 

GG(1,J)=B3NV(I,J) 

GG(l,k)=PART4(I,J) 

CONTINUE 

RETURN 

END 


SUBROUTINE  INVERT  (A,  Y,N,NP) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A (100, 100) ,Y(100,100) ,INDX(100) 
DO  12  1=1, N 
DO  11  J=1,N 
Y(I,J)=0. 

CONTINUE 

Y(I,I)-1. 

CONTINUE 

CALL  LUDCMP(A,N,NP,INDX,D) 

DO  13  J=1,N 

CALL  IUBKSB(A,N,NP,INDX,Y(1,  J) ) 

CONTINUE 

REIURN 

END 


SUBROUTINE  UJDCMP(A,N,NP,INDX,D) 
IMPLICIT  REAL*8(A-H,OZ) 

PARAMETER  (NMAX=100,TTNY=1.0E-20) 
DIMENSION  A(100, 100)  ,INDX(100)  ,W(100) 
D=l. 

DO  12  1=1, N 
AAMAX=0. 


DO  11  J=1,N 

IF  (DABS(A(I,J) )  .GT.AAMAX)  AAMAX=DABS(A(I,J) ) 

11  CONTINUE 

IF  (AAMAX.EQ.O.)  PAUSE  'SINGULAR  MATRIX' 
W(I)=1./AAMAX 

12  CONTINUE 
DO  19  J=1,N 
DO  14  1=1, J-l 
SUM=A(I,  J) 

DO  13  K=l, 1-1 
SUM=SUM-A  ( I ,  K)  *A  (K,  J) 

13  CONTINUE 
A(I, J)=SUM 

14  CONTINUE 
AAMAX=0. 

DO  16  I=J,N 
SUM=A(I,  J) 

DO  15  K=l, J-l 
SUM=SUM-A  ( I ,  K)  *A  (K,  J) 

15  CONTINUE 
A(I,J)=SUM 
DUM=W(I)  *DABS  (SUM) 

IF(DUM.GE.AAMAX)  THEN 
IMAX=I 

AAMAX=DUM 

ENDIF 

16  CONTINUE 
IF(J.NE.IMAX)  THEN 
DO  17  K=1,N 
DUM^A(IMAX,K) 

A(IMAX,K)=A(J,K) 

A(J,K)=DUM 

17  CCNITNUE 
D=-D 

W  (IMAX)  =W  (J) 

ENDIF 

INDX  ( J)  =IMAX 

IF(A(J,J)  .EQ.O. )  A(J, J)=TINY 
IF (J.NE.N)  THEN 
DUM=1./A(J,J) 

DO  18  I=J+1,N 
A(I,  J)=A(I,  J)  *DUM 

18  CONTINUE 
ENDIF 

19  CONTINUE 
RETURN 
END 

C 

C 

C 
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SUBROUTINE  IUBKSB(A,N,NP/INDX,B) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  A(100, 100)  /INDX(100)  ,B(100) 
11=0 

DO  12  1=1 ,N 
LLfINDX(I) 

SUM=B(LL) 

B(LL)=B(I) 

IF  (II.NE.O)  THEN 
DO  11  J=II,I-1 
SUM=SUM-A(I,  J)  *B(J) 

11  CONTINUE 

ELSEIF  (SUM.NE.O.)  THEN 

II=I 

ENDIF 

B(I)=SUM 

12  CONTINUE 

DO  14  I=N, 1,-1 
SUM=B(I) 

IF(I.IT.N)  THEN 
DO  13  J=I+1,N 
SUM=SUM-A(I,  J)  *B(J) 

13  CONTINUE 
ENDIF 

B(I)=SUTVA(I,I) 

14  CONTINUE 
RETURN 
END 
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